Potential Carbon Storage Facilities Near Import/Export Ports

Import and Procedural Functions

import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
import folium
import contextily as cx
import rtree
from zlib import crc32
import hashlib
from shapely.geometry import Point, LineString, Polygon
import numpy as np
from scipy.spatial import cKDTree
from shapely.geometry import Point
from haversine import Unit
from geopy.distance import distance
/Users/jnapolitano/venvs/finance/lib/python3.9/site-packages/geopandas/_compat.py:111: UserWarning: The Shapely GEOS version (3.10.2-CAPI-1.16.0) is incompatible with the GEOS version PyGEOS was compiled with (3.10.1-CAPI-1.16.0). Conversions between both will be slow.
  warnings.warn(

Restrictions

  • Must be near a pipeline terminal

  • Must be Near a petrolium Port

Query Plan

Imports

  • Import LNG terminal Dataa

  • Import well data

Filtering

  • for each well calculate nearest terminal

  • for each well calculate distance from nearest terminal

  • eliminate wells further than 15 miles from a terminal

Data Frame Import

Wells DataFrame

## Importing our DataFrames

gisfilepath = "/Users/jnapolitano/Projects/data/energy/non-active-wells.geojson"


wells_df = gpd.read_file(gisfilepath)

wells_df = wells_df.to_crs(epsg=3857)

Terminal DataFrame

## Importing our DataFrames

gisfilepath = "/Users/jnapolitano/Projects/data/energy/Liquified_Natural_Gas_Import_Exports_and_Terminals.geojson"


terminal_df = gpd.read_file(gisfilepath)

terminal_df = terminal_df.to_crs(epsg=3857)

Eliminating SUSPENDED Terminal from the DataFrame

terminal_df.drop(terminal_df[terminal_df['STATUS'] == 'SUSPENDED'].index, inplace = True)
terminal_df.rename(columns={"NAME": "TERMINAL_NAME"})
terminal_df['TERMINAL_GEO'] = terminal_df['geometry'].copy()
terminal_df.columns
Index(['OBJECTID', 'TERMID', 'NAME', 'ADDRESS', 'CITY', 'STATE', 'ZIP', 'ZIP4',
       'TELEPHONE', 'TYPE', 'STATUS', 'POPULATION', 'COUNTY', 'COUNTYFIPS',
       'COUNTRY', 'LATITUDE', 'LONGITUDE', 'NAICS_CODE', 'NAICS_DESC',
       'SOURCE', 'SOURCEDATE', 'VAL_METHOD', 'VAL_DATE', 'WEBSITE', 'EPA_ID',
       'ALT_NAME', 'OWNER', 'POSREL', 'JRSDCTN', 'CONTYPE', 'IE_PORT',
       'BERTHS', 'STORAGE', 'STORCAP', 'CURRENTCAP', 'APPCAP', 'OPYEAR',
       'IMPEXPCTRY', 'VOLUME', 'PRICE', 'geometry', 'TERMINAL_GEO'],
      dtype='object')

Plotting Terminal by TYPE

terminal_map =terminal_df.explore(
    column="TYPE", # make choropleth based on "PORT_NAME" column
     popup=True, # show all values in popup (on click)
     tiles="Stamen Terrain", # use "CartoDB positron" tiles
     cmap='Set1', # use "Set1" matplotlib colormap
     #style_kwds=dict(color="black"),
     marker_kwds= dict(radius=6),
     #tooltip=['NAICS_DESC','REGION', 'COMMODITY' ],
     legend =True, # use black outline)
     categorical=True,
    )


terminal_map
Make this Notebook Trusted to load map: File -> Trust Notebook

Terminal Impressions

According to the data there is not an export nor import location on The Western Side of the United States.

East Asian import or carbon capture export demands could justfity port development. Another study must be conducted.

Filtering Wells by Terminal Distance in Scipy

Edit

This method does not accuraletly calculate distance. The algorith used below calculates distance on a euclidan plane. In order to calculate a correct answer we must account for sphericity.

I include the code below as reference and a learning opportunity

def ckdnearest(gdA, gdB):

    nA = np.array(list(gdA.geometry.apply(lambda x: (x.x, x.y))))
    nB = np.array(list(gdB.geometry.apply(lambda x: (x.x, x.y))))
    btree = cKDTree(nB)
    dist, idx = btree.query(nA, k=1)
    gdB_nearest = gdB.iloc[idx].drop(columns="geometry").reset_index(drop=True)
    gdf = pd.concat(
        [
            gdA.reset_index(drop=True),
            gdB_nearest,
            pd.Series(dist, name='dist')
        ], 
        axis=1)

    return gdf

c = ckdnearest(wells_df, terminal_df)
c.describe()
index OBJECTID LATITUDE LONGITUDE PERMITNO OPERATORID SURF_LAT SURF_LONG BOT_LAT BOT_LONG ... LATITUDE LONGITUDE BERTHS STORAGE STORCAP CURRENTCAP APPCAP VOLUME PRICE dist
count 1.604460e+05 1.604460e+05 160446.000000 160446.000000 1.604460e+05 1.604460e+05 160446.000000 160446.000000 160432.000000 160446.000000 ... 160446.000000 160446.000000 160446.00000 160446.000000 160446.000000 160446.000000 160446.000000 160446.000000 160446.000000 1.604460e+05
mean 3.596509e+05 3.596519e+05 37.195605 -89.030287 1.980508e+06 6.624020e+07 -63.243732 -176.705395 -949.703169 -955.920831 ... 34.083355 -87.268514 -33.95252 -30.998635 -125.730124 1.641901 -918.544126 -131.529641 -131.638380 7.437940e+05
std 2.806046e+05 2.806046e+05 4.464367 8.648121 3.210840e+06 2.528493e+08 306.472981 269.412702 220.324332 192.546614 ... 4.011683 12.477050 185.97266 186.552966 346.868737 0.518523 272.088477 344.706936 344.554868 6.029022e+05
min 0.000000e+00 1.000000e+00 26.047093 -123.342380 -9.990000e+02 -9.990000e+02 -999.000000 -999.000000 -999.000000 -999.000000 ... 27.889869 -116.847415 -999.00000 -999.000000 -999.000000 0.000000 -999.000000 -999.000000 -999.000000 5.272348e+02
25% 1.392022e+05 1.392032e+05 32.810740 -93.942380 -9.990000e+02 -9.990000e+02 31.861820 -94.036088 -999.000000 -999.000000 ... 30.112960 -93.288224 2.00000 2.000000 7.500000 1.500000 -999.000000 0.000000 0.000000 3.672892e+05
50% 3.128155e+05 3.128165e+05 38.293375 -87.987450 -9.990000e+02 -9.990000e+02 38.067251 -89.171775 -999.000000 -999.000000 ... 31.988522 -88.503111 2.00000 4.000000 10.400000 1.800000 -999.000000 0.000000 0.000000 5.391763e+05
75% 6.171978e+05 6.171988e+05 39.201190 -81.236495 3.305543e+06 -9.990000e+02 39.186490 -81.239420 -999.000000 -999.000000 ... 38.389603 -76.409092 2.00000 7.000000 14.800000 1.800000 -999.000000 9.100000 9.370000 9.642862e+05
max 1.506196e+06 1.506197e+06 50.744220 -75.593800 2.017004e+07 1.044775e+09 50.744220 -75.593800 45.179110 -78.423180 ... 42.392856 -71.058151 2.00000 7.000000 18.000000 4.000000 4.000000 932.000000 11.340000 3.191132e+06

8 rows × 24 columns

nearest_wells_df= wells_df.sjoin_nearest(terminal_df, distance_col="distance_euclidian")
nearest_wells_df.describe()
index OBJECTID_left LATITUDE_left LONGITUDE_left PERMITNO OPERATORID SURF_LAT SURF_LONG BOT_LAT BOT_LONG ... LATITUDE_right LONGITUDE_right BERTHS STORAGE STORCAP CURRENTCAP APPCAP VOLUME PRICE distance_euclidian
count 1.604460e+05 1.604460e+05 160446.000000 160446.000000 1.604460e+05 1.604460e+05 160446.000000 160446.000000 160432.000000 160446.000000 ... 160446.000000 160446.000000 160446.00000 160446.000000 160446.000000 160446.000000 160446.000000 160446.000000 160446.000000 1.604460e+05
mean 3.596509e+05 3.596519e+05 37.195605 -89.030287 1.980508e+06 6.624020e+07 -63.243732 -176.705395 -949.703169 -955.920831 ... 34.083355 -87.268514 -33.95252 -30.998635 -125.730124 1.641901 -918.544126 -131.529641 -131.638380 7.437940e+05
std 2.806046e+05 2.806046e+05 4.464367 8.648121 3.210840e+06 2.528493e+08 306.472981 269.412702 220.324332 192.546614 ... 4.011683 12.477050 185.97266 186.552966 346.868737 0.518523 272.088477 344.706936 344.554868 6.029022e+05
min 0.000000e+00 1.000000e+00 26.047093 -123.342380 -9.990000e+02 -9.990000e+02 -999.000000 -999.000000 -999.000000 -999.000000 ... 27.889869 -116.847415 -999.00000 -999.000000 -999.000000 0.000000 -999.000000 -999.000000 -999.000000 5.272348e+02
25% 1.392022e+05 1.392032e+05 32.810740 -93.942380 -9.990000e+02 -9.990000e+02 31.861820 -94.036088 -999.000000 -999.000000 ... 30.112960 -93.288224 2.00000 2.000000 7.500000 1.500000 -999.000000 0.000000 0.000000 3.672892e+05
50% 3.128155e+05 3.128165e+05 38.293375 -87.987450 -9.990000e+02 -9.990000e+02 38.067251 -89.171775 -999.000000 -999.000000 ... 31.988522 -88.503111 2.00000 4.000000 10.400000 1.800000 -999.000000 0.000000 0.000000 5.391763e+05
75% 6.171978e+05 6.171988e+05 39.201190 -81.236495 3.305543e+06 -9.990000e+02 39.186490 -81.239420 -999.000000 -999.000000 ... 38.389603 -76.409092 2.00000 7.000000 14.800000 1.800000 -999.000000 9.100000 9.370000 9.642862e+05
max 1.506196e+06 1.506197e+06 50.744220 -75.593800 2.017004e+07 1.044775e+09 50.744220 -75.593800 45.179110 -78.423180 ... 42.392856 -71.058151 2.00000 7.000000 18.000000 4.000000 4.000000 932.000000 11.340000 3.191132e+06

8 rows × 25 columns

Calculating Distance in Kilometers from Import/Export Terminal

#df.geopy.distance.distance(coords_1, coords_2).km
#df.apply(lambda row: distance(row['point'], row['point_next']).km if row['point_next'] is not None else float('nan'), axis=1)
# Thanks to https://stackoverflow.com/questions/55909305/using-geopy-in-a-dataframe-to-get-distances

nearest_wells_df['true_distance_km'] = nearest_wells_df.apply(lambda row: distance((row['LATITUDE_left'], row['LONGITUDE_left']), (row['LATITUDE_right'], row['LONGITUDE_right'])).km if row['geometry'] is not None else float('nan'), axis=1)
nearest_wells_df.describe()
index OBJECTID_left LATITUDE_left LONGITUDE_left PERMITNO OPERATORID SURF_LAT SURF_LONG BOT_LAT BOT_LONG ... LONGITUDE_right BERTHS STORAGE STORCAP CURRENTCAP APPCAP VOLUME PRICE distance_euclidian true_distance_km
count 1.604460e+05 1.604460e+05 160446.000000 160446.000000 1.604460e+05 1.604460e+05 160446.000000 160446.000000 160432.000000 160446.000000 ... 160446.000000 160446.00000 160446.000000 160446.000000 160446.000000 160446.000000 160446.000000 160446.000000 1.604460e+05 160446.000000
mean 3.596509e+05 3.596519e+05 37.195605 -89.030287 1.980508e+06 6.624020e+07 -63.243732 -176.705395 -949.703169 -955.920831 ... -87.268514 -33.95252 -30.998635 -125.730124 1.641901 -918.544126 -131.529641 -131.638380 7.437940e+05 592.941775
std 2.806046e+05 2.806046e+05 4.464367 8.648121 3.210840e+06 2.528493e+08 306.472981 269.412702 220.324332 192.546614 ... 12.477050 185.97266 186.552966 346.868737 0.518523 272.088477 344.706936 344.554868 6.029022e+05 467.343457
min 0.000000e+00 1.000000e+00 26.047093 -123.342380 -9.990000e+02 -9.990000e+02 -999.000000 -999.000000 -999.000000 -999.000000 ... -116.847415 -999.00000 -999.000000 -999.000000 0.000000 -999.000000 -999.000000 -999.000000 5.272348e+02 0.454603
25% 1.392022e+05 1.392032e+05 32.810740 -93.942380 -9.990000e+02 -9.990000e+02 31.861820 -94.036088 -999.000000 -999.000000 ... -93.288224 2.00000 2.000000 7.500000 1.500000 -999.000000 0.000000 0.000000 3.672892e+05 311.703128
50% 3.128155e+05 3.128165e+05 38.293375 -87.987450 -9.990000e+02 -9.990000e+02 38.067251 -89.171775 -999.000000 -999.000000 ... -88.503111 2.00000 4.000000 10.400000 1.800000 -999.000000 0.000000 0.000000 5.391763e+05 421.824879
75% 6.171978e+05 6.171988e+05 39.201190 -81.236495 3.305543e+06 -9.990000e+02 39.186490 -81.239420 -999.000000 -999.000000 ... -76.409092 2.00000 7.000000 14.800000 1.800000 -999.000000 9.100000 9.370000 9.642862e+05 799.595274
max 1.506196e+06 1.506197e+06 50.744220 -75.593800 2.017004e+07 1.044775e+09 50.744220 -75.593800 45.179110 -78.423180 ... -71.058151 2.00000 7.000000 18.000000 4.000000 4.000000 932.000000 11.340000 3.191132e+06 2390.537825

8 rows × 26 columns

Filtering Wells within 50 KM of a Terminal

filtered_wells = nearest_wells_df.loc[nearest_wells_df['true_distance_km'] < 50].copy()
filtered_wells.describe()
index OBJECTID_left LATITUDE_left LONGITUDE_left PERMITNO OPERATORID SURF_LAT SURF_LONG BOT_LAT BOT_LONG ... LONGITUDE_right BERTHS STORAGE STORCAP CURRENTCAP APPCAP VOLUME PRICE distance_euclidian true_distance_km
count 2.654000e+03 2.654000e+03 2654.000000 2654.000000 2654.000000 2654.0 2654.000000 2654.000000 2652.000000 2654.000000 ... 2654.000000 2654.000000 2654.000000 2654.000000 2654.000000 2654.000000 2654.000000 2654.000000 2654.000000 2654.000000
mean 1.554318e+05 1.554328e+05 30.005777 -93.585830 16514.931047 -999.0 29.230172 -94.272273 -908.678451 -919.670737 ... -93.524363 -375.983798 -375.082894 -370.543293 0.973576 -350.589940 -374.631349 -376.703112 24095.468770 20.848131
std 3.226352e+05 3.226352e+05 0.377108 0.665728 90735.288423 0.0 28.244840 24.858548 291.080631 255.761340 ... 0.630421 485.299749 486.001893 489.542071 1.090584 478.695054 488.693585 484.744059 14264.583731 12.326692
min 2.500000e+01 2.600000e+01 27.616626 -97.747110 -999.000000 -999.0 -999.000000 -999.000000 -999.000000 -999.000000 ... -97.274403 -999.000000 -999.000000 -999.000000 0.000000 -999.000000 -999.000000 -999.000000 527.234827 0.454603
25% 1.038925e+04 1.039025e+04 29.986090 -93.621635 -999.000000 -999.0 29.985485 -93.621915 -999.000000 -999.000000 ... -93.337457 -999.000000 -999.000000 -999.000000 0.000000 -999.000000 -999.000000 -999.000000 12302.653045 10.625148
50% 5.860500e+04 5.860600e+04 30.032030 -93.441460 -999.000000 -999.0 30.031910 -93.442460 -999.000000 -999.000000 ... -93.336491 2.000000 3.000000 9.000000 0.710000 1.410000 0.000000 0.000000 26936.447255 23.256971
75% 1.096650e+05 1.096660e+05 30.251400 -93.338690 -999.000000 -999.0 30.251395 -93.338887 -999.000000 -999.000000 ... -93.329114 2.000000 3.000000 11.000000 1.800000 4.000000 0.000000 0.000000 32607.124822 28.267032
max 1.503142e+06 1.503143e+06 30.561100 -88.052220 833497.000000 -999.0 30.561100 -92.782800 30.077654 -93.845118 ... -88.503111 2.000000 5.000000 18.000000 4.000000 4.000000 932.000000 11.340000 58104.197789 49.997185

8 rows × 26 columns

Map of Wells within 50 km of an Import/Export Terminal by Type

filtered_wells.explore(
    column="STATUS_left", # make choropleth based on "PORT_NAME" column
     popup=True, # show all values in popup (on click)
     tiles="Stamen Terrain", # use "CartoDB positron" tiles
     cmap='Set1', # use "Set1" matplotlib colormap
     #style_kwds=dict(color="black"),
     marker_kwds= dict(radius=6),
     #tooltip=['NAICS_DESC','REGION', 'COMMODITY' ],
     legend =True, # use black outline)
     categorical=True,)
Make this Notebook Trusted to load map: File -> Trust Notebook